function[den_grid_tot,vp_grid_tot,vs_grid_tot, corr_cell_tot]=...
    select_best_AVA_regression(angles, rangeOfAngles, weights_best,corr_cell_tot,layers,den_grid_tot,vp_grid_tot,vs_grid_tot)
%% EVALUATES THE CORRELATION FOR EACH TRACE AS A FUNCTION OF THE ANGLE
% CREATES BEST FOR EACH ITERATION AND AND GLOBAL BEST FOR noIt > 1
% GLOBAL FROM noIt-1 IS INCLUDED IN THE LS EVALUATION
%
% INPUT TYPE - ALL REALIZATION FROM noIt ITERATION
%
% corr_cell [Nx1] - CORRELATION CELL - only the ones that correspond to
%                   noIt iteration
% angles - ANGLES THAT BELONG TO THE GATHER
% noLayers - NUMBER OF LAYERS USED FOR THE CORR COEF CALCULATION
% sy_cell - SEISMIC CELL WITH SYNTHETIC FROM noIt ITERATION
% D_grid - DENSITY SIMULATIONS FOR THE noIt INTERATION
% Vp_grid - Vp SIMULATIONS FOR THE noIt INTERATION
% Vs_grid - Vs SIMULATIONS FOR THE noIt INTERATION
% Rzero_grid - GRID WITH Rzero RESULTS FOR THE ENTIRE ITERATION
% G_grid - GRID WITH G RESULTS FOR THE ENTIRE ITERATION
%
% best_cc_it [cell] - GATHER FROM corr_cell WITH THE BEST CC. ONE PER
%                         ITERATION
% best_W [noLayers x 1] - REALIZATION NUMBER WITH THE BEST LAYER
% W [noLayers x noSim] - WEIGHTS. ROWS - LAYER NUMBER; COL. - REALIZATION NUMBER
% interp - VALUES FOR WEIGTHED SELECTION:
%   - interp(N,1) - WEIGHTS FOR SLOPE
%   - interp(N,2) - WEIGHTS FOR C
%   - interp(N,3) - WEIGHTS FOR LS
% y_int - INTERPOLATED Y VALUES
% p [N x 2] - RETURN THE COEFFICIENTS OF THE LINEAR FIT
% LS [N x 1]-  LEAST SQUARES FOR EACH SIM.
%
% FUNCTIONS
% grid2gslib.m
%
% MODIFICATIONS
% 12 DEZ 2013 - REMOVE CELL FOR VARIABLES p AND W
% 27 NOV 2013 - TEST FOR 1 ANGLE ONLY - BEHAVE LIKE GSI
% 12 JUL 2013 - BUG FIXED - INTERSECT VALUE WAS NOT WELL INTERP.
% 29 JUN 2013 - CHANGES TO STORE ONLY CURRENT SIM AND BEST
% 12 MAR 2013 - Rzero AND G OUTPUTS WERE CHANGED IN ORDER
% 08 OCT 2012 - NUMERIAL APPROXIMATION ERROR (SINGLE) CORRECTED
% 04 OCT 2012 - NUMERIAL APPROXIMATION ERROR (SINGLE) CORRECTED
% 28 SEP 2012 - NEW CODE FOR CHECK INCREASING ORDER BEFORE INTERPOLATION
% 29 AUG 2012 - PERFORMANCE UPDATE - p & LS
% 23 AUG 2012 - polyfit FUNCTION REMOVED FOR \ - PERFORMANCE INCREASE
% 03 AUG 2012 - SAVES best_vp & best_vs; REORDER SIMULATIONS PER COL.
%               CHANGED
% 31 JUL 2012 - PERFORMANCE UPDATE
% 27 APR 2012 - FIXED FOR VARIABLE SIZE LAYERS
% 24 APR 2012 - BEST GLOBAL FROM IT-1 NOT CONSIDERED
% 22 APR 2012 - LAYERS OF RANDOM SIZE
%
% LA - JUL 2019
%

%% INITIALIZE AND ALLOCATE

size_grid            = size(den_grid_tot{1,1});
corr_cell_tot{1,1}   = cell2mat(corr_cell_tot{1,1});
corr_cell_tot{2,1}   = cell2mat(corr_cell_tot{2,1});
bestCC               = zeros(size(corr_cell_tot{1,1},1),size(corr_cell_tot{1,1},2));

interp = cell(1,4);
for n = 1:size(interp,2)
    interp{1,n} = zeros(size(corr_cell_tot{1,1},1),2);
end

sizeAngles = size(angles,2);
sizeRangeofAngles = size(rangeOfAngles(1,1):rangeOfAngles(1, end));
%% MAIN
for l = 1:size_grid(2)
          
    A_corr      = [ones(sizeRangeofAngles(1,2),1), angles(sizeRangeofAngles(1,1):sizeRangeofAngles(1,2))']; % CREATE ARRAY FOR POLYNOMIAL FIT
    
    pCurrent    = fliplr((A_corr\permute(corr_cell_tot{1,1}(:,((l-1) * sizeAngles) + 1:(l-1) * sizeAngles + sizeRangeofAngles(1,2)),[2 1]))'); % COMPUTES THE MEAN p(1) AND THE INTESECT p(2) FOR THE CORRELATION VALUES
    pBest       = fliplr((A_corr\permute(corr_cell_tot{2,1}(:,((l-1) * sizeAngles) + 1:(l-1) * sizeAngles + sizeRangeofAngles(1,2)),[2 1]))'); % COMPUTES THE MEAN p(1) AND THE INTESECT p(2) FOR THE CORRELATION VALUES
    
%FT     yIntCurrent = pCurrent(:,1) .* repmat(angles(sizeRangeofAngles(1,1):sizeRangeofAngles(1,2)), size(corr_cell_tot{1,1},1),1) + pCurrent(:,2);
%FT     yIntBest    = pBest(:,1)    .* repmat(angles(sizeRangeofAngles(1,1):sizeRangeofAngles(1,2)), size(corr_cell_tot{2,1},1),1) + pBest(:,2);
    
    % compute LS
%FT     interp{1,3}(:,1) = sum((corr_cell_tot{1,1}(:,((l-1) * sizeAngles) + 1:(l-1) * sizeAngles + sizeRangeofAngles(1,2)) - yIntCurrent).^2,2);
%FT     interp{1,3}(:,2) = sum((corr_cell_tot{2,1}(:,((l-1) * sizeAngles) + 1:(l-1) * sizeAngles + sizeRangeofAngles(1,2)) - yIntBest).^2,2);
    
    % CHECK ZERO FOR SLOPE
    interp{1,1}(pCurrent(:,1) == 0,1)   = 1;
    interp{1,1}(pCurrent(:,1) ~= 0,1)   = pCurrent(pCurrent(:,1) ~= 0,1);
    interp{1,1}(pBest(:,1)    == 0,2)   = 1;
    interp{1,1}(pBest(:,1)    ~= 0,2)   = pBest(pBest(:,1) ~= 0,1);

    % CHECK ZERO FOR C
    interp{1,2}(pCurrent(:,2) == 0,1)   = 0.00001;
    interp{1,2}(pCurrent(:,2) ~= 0,1)   = pCurrent(pCurrent(:,2) ~= 0, 2);
    interp{1,2}(pBest(:,2) == 0,2)      = 0.00001;
    interp{1,2}(pBest(:,2) ~= 0,2)      = pBest(pBest(:,2) ~= 0, 2);
    
    % CHECK ZERO FOR LS
    interp{1,3}(interp{1,3}(:,1) == 0, 1) = 1;
    interp{1,3}(interp{1,3}(:,2) == 0, 2) = 1;
        
    % CHECK ZERO FOR MEAN(CC)
    temp    = mean(corr_cell_tot{1,1}(:,((l-1) * sizeAngles) + 1:(l-1) * sizeAngles + sizeRangeofAngles(1,2)), 2);
    interp{1,4}(temp == 0, 1) = 0.00001;
    interp{1,4}(temp ~= 0, 1) = temp(temp ~= 0);
    
    temp    = mean(corr_cell_tot{2,1}(:,((l-1) * sizeAngles) + 1:(l-1) * sizeAngles + sizeRangeofAngles(1,2)), 2);
    interp{1,4}(temp == 0, 2) = 0.00001;
    interp{1,4}(temp ~= 0, 2) = temp(temp ~= 0);

    % W is the weighted sum of each parameter in the interpolation 
    W = weights_best(1,1).*(1./abs(interp{1,1})./repmat(max((1./abs(interp{1,1})),[],2),1,2))+...
        weights_best(1,2).*(interp{1,2}./repmat(max(abs(interp{1,2}),[],2),1,2))+...
        weights_best(1,3).*(1./interp{1,3}./repmat(max(abs(1./interp{1,3}),[],2),1,2))+...
        weights_best(1,4).*(interp{1,4}./repmat(max(abs(interp{1,4}),[],2),1,2));
    
    % find best positions
    [~,best_W] = max(W,[],2); %returns max for each col (simulation)

    % build BCC gathers
    bestCC(best_W(:,1) == 1,((l-1) * sizeAngles) + 1:(l-1) * sizeAngles + sizeAngles) = corr_cell_tot{1,1}(best_W(:,1) == 1,((l-1) * sizeAngles) + 1:(l-1) * sizeAngles + sizeAngles);
    bestCC(best_W(:,1) == 2,((l-1) * sizeAngles) + 1:(l-1) * sizeAngles + sizeAngles) = corr_cell_tot{2,1}(best_W(:,1) == 2,((l-1) * sizeAngles) + 1:(l-1) * sizeAngles + sizeAngles);
 
    indLayer        = 1;
    bestWReshape    = reshape(best_W,size(layers,1),size_grid(1));
    
    for n = 1: size(layers,1)
        
        % best volumes for elastic properties
        den_grid_tot{2,1}(bestWReshape(n,:)==1, l, indLayer:indLayer+layers(n)-1) = den_grid_tot{1,1}(bestWReshape(n,:)==1, l, indLayer:indLayer+layers(n)-1);
        den_grid_tot{2,1}(bestWReshape(n,:)==2, l, indLayer:indLayer+layers(n)-1) = den_grid_tot{2,1}(bestWReshape(n,:)==2, l, indLayer:indLayer+layers(n)-1);

        vp_grid_tot{2,1}(bestWReshape(n,:)==1, l, indLayer:indLayer+layers(n)-1)  = vp_grid_tot{1,1}(bestWReshape(n,:)==1,  l, indLayer:indLayer+layers(n)-1);
        vp_grid_tot{2,1}(bestWReshape(n,:)==2, l, indLayer:indLayer+layers(n)-1)  = vp_grid_tot{2,1}(bestWReshape(n,:)==2,  l, indLayer:indLayer+layers(n)-1);
        
        vs_grid_tot{2,1}(bestWReshape(n,:)==1, l, indLayer:indLayer+layers(n)-1)   = vs_grid_tot{1,1}(bestWReshape(n,:)==1,  l, indLayer:indLayer+layers(n)-1);
        vs_grid_tot{2,1}(bestWReshape(n,:)==2, l, indLayer:indLayer+layers(n)-1)   = vs_grid_tot{2,1}(bestWReshape(n,:)==2,  l, indLayer:indLayer+layers(n)-1);
  
        indLayer    = indLayer+layers(n,1);
    end

end

corr_cell_tot{1,1} = mat2cell(corr_cell_tot{1,1}, repmat(size(layers,1),size_grid(1),1),  repmat(sizeAngles,size_grid(2),1));
corr_cell_tot{2,1} = mat2cell(bestCC, repmat(size(layers,1),size_grid(1),1),  repmat(sizeAngles,size_grid(2),1));        
    
end




